%% term_grad_P_grad_Q
function matrix_A = term_grad_P_grad_Q(mesh, fespace)
Nb = fespace.Nb;
Nlb = fespace.Nlb;
rows = zeros(mesh.Ne*Nlb*Nlb,1);
cols = rows;
val = rows;
flag = 1;
for e = 1:mesh.Ne
    V = mesh.node(mesh.elem(e,:),:);
    [int_points_loc, int_weights_loc] = integral_loc(V, "3d");
    phix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dx");
    phiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dy");
    phiz = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dz");
    for i = 1:Nlb
        for j = 1:Nlb
            f = phix(:,j).*phix(:,i) + phiy(:,j).*phiy(:,i) + phiz(:,j).*phiz(:,i);
            r = int_weights_loc * f;
            rows(flag) = fespace.elem(e,i);
            cols(flag) = fespace.elem(e,j);
            val(flag) = r;
            flag = flag + 1;
        end
    end
end
matrix_A = sparse(rows, cols, val, Nb, Nb);
end